El avance en las tecnologías de secuenciación en los últimos años, ha permitido estudiar mayor número de secuencias genéticas provenientes de distintos organismos. Asímismo, los avances en el almacenamiento de la información secuenciada ha favorecido el análisis masivo de datos.
El sufijo “omica” refiere el estudio completo de dicho nivel molecular.
La secuenciación de transcriptomas completos permite conocer los genes que expresa una célula o un conjunto de células en un momento determinado y bajo condiciones particulares.
Algunas de las preguntas que pueden responderse con el análisis de transcriptomas son:
La técnica que permite realizar la secuenciación de los transcriptomas es la secuenciación del RNA (RNA-seq).
Una vez secuenciadas las bibliotecas se genera una cantidad masiva de datos ¿qué a hacer con toda la información obtenida?
La cantidad de datos que se generan a partir de la secuenciación de RNA es inmensa. El análisis de expresión diferencial requiere seguir el siguiente pipeline o flujo de trabajo:
pipeline
Al finalizar el proceso de estimación de la abundancia, RSEM nos devuelve los siguientes archivos por muestra:
## t05.genes.results
## t05.isoforms.results
## t05.txt
Y tienen el siguiente contenido de archivos de texto plano:
## gene_id transcript_id(s) length effective_length expected_count TPM FPKM
## ENSG00000000003 ENST00000373020,ENST00000494424,ENST00000496771,ENST00000612152,ENST00000614008 2156.34 1984.51 756.00 33.08 25.14
## ENSG00000000005 ENST00000373031,ENST00000485971 940.50 768.71 0.00 0.00 0.00
## ENSG00000000419 ENST00000371582,ENST00000371584,ENST00000371588,ENST00000413082,ENST00000466152,ENST00000494752 1077.25 905.41 517.00 49.58 37.68
Se requieren importar los datos generados (genes.results o .transcripts.results) por RSEM a R. Para ello se requiere de la librería de tximport [(??? ctbTximport2017). Esta librería importa archivos de cuentas generados por:
Se requiere de una tabla de metadatos:
## Sample Cell_Line Condition Replicate Name
## 1 t03 A CT 1 A_CT_1
## 2 t05 A Venom 2 A_Venom_2
## 3 t07 A CT 2 A_CT_2
## 4 t013 A Venom 1 A_Venom_1
## 5 t026 A Venom 3 A_Venom_3
## 6 t028 A CT 3 A_CT_3
Y crear una ruta hacia la ubicación de los archivos:
##Se requiere crear una ruta al directorio en donde se encuentran los archivos
files <- file.path("../rsem/rsem/", samples$Sample, paste0(samples$Sample, ".genes.results"))
names(files) <- samples$Name
##Se corrobora que la ruta fue creada y contiene los nombres de los archivos correctamente
data.frame(files = files, Exists = file.exists(files))## files Exists
## A_CT_1 ../rsem/rsem//t03/t03.genes.results TRUE
## A_Venom_2 ../rsem/rsem//t05/t05.genes.results TRUE
## A_CT_2 ../rsem/rsem//t07/t07.genes.results TRUE
## A_Venom_1 ../rsem/rsem//t013/t013.genes.results TRUE
## A_Venom_3 ../rsem/rsem//t026/t026.genes.results TRUE
## A_CT_3 ../rsem/rsem//t028/t028.genes.results TRUE
Cuidado: Revisar bien el orden (jerarquía) y nombre de los directorios y archivos.
Se importan los archivos empleando tximport:
## reading in files with read_tsv
## 1 2 3 4 5 6
## List of 4
## $ abundance : num [1:58735, 1:6] 24.92 0 72.28 3.18 15.73 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:58735] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
## .. ..$ : chr [1:6] "A_CT_1" "A_Venom_2" "A_CT_2" "A_Venom_1" ...
## $ counts : num [1:58735, 1:6] 998 0 1274 199 715 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:58735] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
## .. ..$ : chr [1:6] "A_CT_1" "A_Venom_2" "A_CT_2" "A_Venom_1" ...
## $ length : num [1:58735, 1:6] 2051 760 903 3206 2328 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:58735] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
## .. ..$ : chr [1:6] "A_CT_1" "A_Venom_2" "A_CT_2" "A_Venom_1" ...
## $ countsFromAbundance: chr "no"
## [1] "matrix"
## A_CT_1 A_Venom_2 A_CT_2 A_Venom_1 A_Venom_3 A_CT_3
## ENSG00000000003 998.00 756.0 1237.00 970.00 1663.00 1279.0
## ENSG00000000005 0.00 0.0 0.00 0.00 0.00 0.0
## ENSG00000000419 1274.00 517.0 1810.00 610.00 1165.00 1549.0
## ENSG00000000457 199.05 177.7 278.53 140.14 173.92 369.4
## ENSG00000000460 714.95 383.3 784.47 611.86 1494.08 796.6
## ENSG00000000938 0.00 0.0 0.00 0.00 0.00 2.0
Para el análisis de expresión diferencial requerimos que los valores de las cuentas sean números enteros. Recodificar los valores a integer y convertir la matriz de cuentas a un data frame:
storage.mode(counts) <- "integer"
##Recodificar la matriz de cuentas a data frame
counts <- as.data.frame(counts)
class(counts)## [1] "data.frame"
## [1] "A_CT_1" "A_Venom_2" "A_CT_2" "A_Venom_1" "A_Venom_3" "A_CT_3"
Filtrado de cuentas Filtrar los datos es muy útil. Los genes cuyos valores de cuentas son igual a 0 no aportan ninguna información interesante al análisis. Además, eliminar genes con bajos valores de cuentas nos evita obtener resultados falsos positivos.
##Seleccionando genes con al menos 5 cuentas por millón (cpm) en al menos 2 muestras
keep <- rowSums(cpm(counts) >= 5) >=2
table(keep)## keep
## FALSE TRUE
## 46685 12050
Es recomendable no realizar un filtrado astringente, de lo contrario se recuperan solamente genes “housekeeping”.
Tanto para edgeR como DESeq2 requerimos almacenar las cuentas y algunos metadatos en un objeto de forma de lista.
Para construir la lista de edgeR requerimos indicar cuáles son los grupos que se compararán y el número de réplicas por grupo:
## [1] "A_CT_1" "A_CT_2" "A_CT_3" "A_Venom_1" "A_Venom_2" "A_Venom_3"
##Usar los nombres de las columnas de la matriz de cuentas para seleccionar los dos grupos experimentales
groups <- factor(sub("..$", "", colnames(counts)))
table(groups)## groups
## A_CT A_Venom
## 3 3
Crear la lista empleando la matriz de cuentas y los grupos:
## Formal class 'DGEList' [package "edgeR"] with 1 slot
## ..@ .Data:List of 3
## .. ..$ : int [1:12050, 1:6] 998 1274 199 714 0 4761 984 1218 313 1315 ...
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:12050] "ENSG00000000003" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" ...
## .. .. .. ..$ : chr [1:6] "A_CT_1" "A_CT_2" "A_CT_3" "A_Venom_1" ...
## .. ..$ :'data.frame': 6 obs. of 3 variables:
## .. .. ..$ group : Factor w/ 2 levels "A_CT","A_Venom": 1 1 1 2 2 2
## .. .. ..$ lib.size : num [1:6] 25791745 35428109 38155734 17030975 15049033 ...
## .. .. ..$ norm.factors: num [1:6] 1 1 1 1 1 1
## .. ..$ :'data.frame': 12050 obs. of 1 variable:
## .. .. ..$ genes: chr [1:12050] "ENSG00000000003" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" ...
Posteriormente los datos son normalizados empleando el método de los TMM (weighted Trimmed Mean of M-values). Para ello edgeR busca y elimina los valores atípicos (expresión absoluta de muestra como expresión relativa entre muestra) y calcula los factores de normalización. En este artículo encontrarán de manera detallada la explicación de la normalización por TMM (Evans, Hardin, and Stoebel 2018).
## group lib.size norm.factors
## A_CT_1 A_CT 25791745 1.1130004
## A_CT_2 A_CT 35428109 1.0882614
## A_CT_3 A_CT 38155734 1.0922808
## A_Venom_1 A_Venom 17030975 0.9605541
## A_Venom_2 A_Venom 15049033 0.9113970
## A_Venom_3 A_Venom 36868378 0.8633914
Es importante inspeccionar los datos y verificar que estén correctamente normalizados. Para ello, se grafica la expresión absoluta vs expresión relativa para cada muestra:
La consistencia de las réplicas la podemos verificar mediante un análisis de componentes principales (PCA) o calculando la correlación (Pearson) que existe entre las muestras:
p <- cpm(edgeRlist$counts, log = T)
p <- pca(p)
biplot(p, lab = colnames(edgeRlist$counts), pointSize = 10,
title = "PCA")cormat <- cor(cpm(edgeRlist$counts, log = T))
pheatmap(cormat, border_color = NA, main = "Correlation of replicates")Para el análisis de expresión diferencial requerimos:
En esta matriz le vamos a indicar a edgeR cuáles muestras corresponden a un grupo o condición experimental. Dado que algunos experimentos pueden tener diseños muy complejos (una muestra pertenece a un tipo celular y a un tratamiento) empleamos la función interna de R, modelado de matrices:
design <- model.matrix(~0+edgeRlist$samples$group)
##El término ~0 le indica a la función no incluir una columna de intersecciones y solamente incluir tantas columnas como grupos en nuestro diseño experimental
colnames(design) <- levels(edgeRlist$samples$group)
design## A_CT A_Venom
## 1 1 0
## 2 1 0
## 3 1 0
## 4 0 1
## 5 0 1
## 6 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`edgeRlist$samples$group`
## [1] "contr.treatment"
Como edgeR ajusta los datos a un modelo de distribución bi-nomial negativo (parecido a Poisson) se requiere calcular un parámetro adicional de dispersión de los datos.
edgeR calcula la dispersión a tres niveles:
edgeR en los análisis de expresión diferencial puede obviar las comparaciones entre los grupos experimentales. Sin embargo, es bueno contar con una matriz de contrastes para indicarle a edgeR cuáles son las comparaciones que queremos hacer entre los datos.
## Contrasts
## Levels Cells
## A_CT -1
## A_Venom 1
En este caso, solo tenemos dos condiciones experimentales “Células A” y “Células B”. Con esta matriz de contrastes indicamos que la comparación será buscando los genes diferencialmente expresados en B con respecto a A. Sin embargo, si tuvieramos mayor cantidad de grupos experimentales los podemos comparar, por ejemplo:
Células CT = X, Y, Z
Células Tx = P, Q, R
makeContrast("Cells" = "(P + Q + R)/3 - (X, Y, Z)/3")
El objetivo es que los coeficientes sumen 0
Los datos deben ajustarse a un modelo lineal bi-nomial negativo. Para ello se utilizará la función glmQLfit con la cual se tiene un control más robusto del “error rate”.
Se realiza la comparación para obtener los genes diferencialmente expresados entre una condición y otra. Con la función glmQLFTest, estamos probando la hipótesis nula “el valor de |lfc| del gen A es igual a 0”. Por lo tanto, los valores de p y p.adj calculados están hechos con respecto a dicha hipótesis nula.
qlf.BvsA <- glmQLFTest(fit, ##Objeto en forma de lista con los datos ajustados a un modelo bi-nomial negativo
contrast = contrast[, "Cells"]) ##Matriz de contrastes
##Obtenemos los DEG con expresión distinta de 1, valor de p menor a 0.05, corrigiendo el valor de p por el método de Benjamini-Hochberg
deg.BvsA <- decideTestsDGE(qlf.BvsA, p.value = 0.05, adjust.method = "BH", lfc = 0)
table(deg.BvsA)## deg.BvsA
## -1 0 1
## 3307 5244 3499
##Guardamos los resultados en un data.frame
DEG.BvsA <- DEGResults(qlf.BvsA)
##Añadimos una columna extra a los resultados anteriores indicando la condición de expresión respecto al lfc y al FDR
DEG.BvsA <- edgeResults(DEG.BvsA, logfc = 0, padj = 0.05)Los datos los visualizamos graficando un “Volcano plot”
Un procedimiento muy común es seleccionar aquellos genes cuya expresión difeerencial haya sido significativa y que cumplan con un valor de lfc. Por ejemplo: |lfc| > 1 (es decir genes con cuya expresión es > 2 o < 0.5 respecto al CT)
significant.1 <- DEG.BvsA %>% filter(logFC > 1 & FDR < 0.05 | logFC < -1 & FDR < 0.05)
paste(length(significant.1$genes), "is the number of significant genes")## [1] "4045 is the number of significant genes"
¡Cuidado! Filtrar los resultados de esta manera es incorrecto. Recordar que en el análisis de expresión diferencial probamos la hipótesis nula “El |lfc| del gen A es igual 0” y los valores de p y FDR corresponden a dicha prueba. El asumir que estos genes filtrados tienen |lfc| > 1 con el valor de FDR previamente calculado provoca un sesgo. Favorece genes con baja expresión y alta variabilidad, invalidando el poder estadístico de la prueba. En otras palabras, incrementamos el riesgo de captar resultados falsos positivos. Aquí encontrarán un artículo en donde se explica detalladamente las implicaciones de realizar este tipo erroneo de filtrados (Love, Huber, and Anders 2014).
Para resolver este problema, tanto edgeR como DESeq2 implementan funciones en donde se prueba la hipótesis nula “el |lfc| del gen A es distinto a x”:
qlf.BvsA.lfc1 <- glmTreat(fit,
contrast = contrast[, "Cells"],
lfc = 1)
deg.BvsA.lfc1 <- decideTestsDGE(qlf.BvsA.lfc1, p.value = 0.05, adjust.method = "BH", lfc = 1)
table(deg.BvsA.lfc1)## deg.BvsA.lfc1
## -1 0 1
## 1198 9405 1447
De esta manera podemos seleccionar los genes que estadísticamente tienen |lfc| > 1 y robustecer nuestros resultados.
DEG.BvsA.lfc1 <- DEGResults(qlf.BvsA.lfc1)
DEG.BvsA.lfc1 <- edgeResults(DEG.BvsA.lfc1, logfc = 1, padj = 0.05)Visualizamos estos datos en un nuevo volcano plot
significant.genes <- DEG.BvsA.lfc1 %>% filter(FDR < 0.05 & logFC > 1 | FDR < 0.05 & logFC < -1)
paste("The number of significant genes with |lfc| > 1 is", length(significant.genes$genes))## [1] "The number of significant genes with |lfc| > 1 is 2645"
Finalmente, los los genes diferencialmente expresados podemos emplearlos para crear mapas de calor “heatmaps”:
##Obtener los nombres o ids de los genes con expresión diferencial
significant.ids <- significant.genes$genes
##Crear una matriz de cuentas normalizadas por cpm empleando las cuentas que se encuentran guardadas en el objeto edgeRlist
significant.cpm <- cpm(edgeRlist$counts, log = T)
##Cortamos los genes con expresión significativa
significant.cpm <- significant.cpm[significant.ids, ]
##Generamos una paleta de colores
RedBlackGreen <- maPalette(low = "green", high = "red", mid = "black")
##Obtenemos el heatmap
pheatmap(significant.cpm,
border_color = NA,
color = RedBlackGreen,
show_rownames = F,
scale = "row",
angle_col = 0)Una vez obtenida la lista de genes con expresión diferencial significativa se procede a realizar análisis de representación de vías (ORA o GSEA) para darle una explicación biloógica a los resultados obtenidos.
En DESeq2 requerímos crear una tabla de metadatos (distinta a la que empleada para importar los datos) con la siguiente información:
##Crear vectores indicando el nombre de las muestras (líneas celulares) y el número de réplicas en cada condición
cell_lines <- c(rep("A", 6))
##Crear un vector indicando las dos condiciones experimentales a comparar
condition <- c(rep("control", 3), rep("tx", 3))
##En un df alojar los datos de las líneas celulares y las condiciones experimentales a las que pertenecen
meta_cells <- data.frame(cell_lines, condition)
##Especificar el nombre de las filas de acuerdo al nombre de las columnas de la matriz de cuentas
rownames(meta_cells) <- names(counts_DESeq)
meta_cells## cell_lines condition
## A_CT_1 A control
## A_CT_2 A control
## A_CT_3 A control
## A_Venom_1 A tx
## A_Venom_2 A tx
## A_Venom_3 A tx
Al igual que en edgeR, DESeq2 guarda o aloja los datos de las cuentas y metadatos en un objeto con clase de lista:
DDS_list <- DESeqDataSetFromMatrix(countData = counts_DESeq, ## Datos de las cuentas crudas
colData = meta_cells, ## Tabla de metadatos
design = ~condition) ## Grupos experimentales a los que pertence cada muestraEn este punto tenemos que especificarle a DESeq cuál es el grupo de referencia para realizar las comparaciones:
Para filtrar los genes con expresión baja, DESeq emplea una función interna de normalización para que los datos de las cuentas sean comparables entre ellos:
##Usamos los mismos criterios que en edgeR
filtered <- rowSums(counts(DDS_list) >= 5) >= 2
table(filtered)## filtered
## FALSE TRUE
## 39421 19314
Las muestras son normalizadas empleando la función DESeq. Esta función realiza tres pasos en un solo comando:
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
Finalmente, obtenemos los genes diferencialmente expresados. En este caso vamos a probar la hipótesis nula “El lfc del gen A es igual a 0”:
##
## out of 19314 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 4932, 26%
## LFC < 0 (down) : 5068, 26%
## outliers [1] : 20, 0.1%
## low counts [2] : 0, 0%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
Visualizamos los datos en un volcano plot:
##convertimos los resultados a un df
DEG.VenomvsCT <- as.data.frame(DEG.VenomvsCT) %>%
rownames_to_column(var = "ensgene_id")
DEG.VenomvsCT <- DESEqResults(DEG.VenomvsCT, logfc = 0, FDR = 0.05)
volcano_DESeq2(DEG.VenomvsCT, lfc = 0, FDR = 0.05)Y de manera similar a edgeR, generamos un heatmap con los genes diferencialmente expresados:
significant.genes.2 <- DEG.VenomvsCT %>% filter(log2FoldChange > 0 & padj < 0.05 | log2FoldChange < 0 & padj < 0.05)
significant.ids.2 <- significant.genes.2$ensgene_id
significant.counts <- counts(DDS_list, normalized = T)
significant.counts <- significant.counts[significant.ids.2, ]
pheatmap(significant.counts,
border_color = NA,
color = RedBlackGreen,
show_rownames = F,
scale = "row",
angle_col = 0)La lista de genes diferencialmente expresados obtenida por edgeR o DESeq2 no aportan información sobre los procesos biológicos que pueden verse afectados en las células tratadas comparado con el control. Para ello, se requiere realizar análisis que le den un contexto biológico a las listas obtenidas. Por ejemplo los análisis de vías son muy útiles para contextualizar los resultados:
En este artículo encontrarán la información detallada de los análisis de vías, sus fortalezas, debilidades y alcances de cada método (García-Campos, Espinal-Enríquez, and Hernández-Lemus 2015).
Dobin, Alexander, Carrie A. Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R. Gingeras. 2013. “STAR: Ultrafast Universal RNA-Seq Aligner.” Bioinformatics 29 (1): 15–21. https://doi.org/10.1093/bioinformatics/bts635.
Evans, Ciaran, Johanna Hardin, and Daniel M. Stoebel. 2018. “Selecting Between-Sample RNA-Seq Normalization Methods from the Perspective of Their Assumptions.” Briefings in Bioinformatics 19 (5): 776–92. https://doi.org/10.1093/bib/bbx008.
García-Campos, Miguel A., Jesús Espinal-Enríquez, and Enrique Hernández-Lemus. 2015. “Pathway Analysis: State of the Art.” Frontiers in Physiology 6: 383. https://doi.org/10.3389/fphys.2015.00383.
Li, Bo, and Colin N. Dewey. 2011. “RSEM: Accurate Transcript Quantification from RNA-Seq Data with or Without a Reference Genome.” BMC Bioinformatics 12 (August): 323. https://doi.org/10.1186/1471-2105-12-323.
Love, Michael I, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” Genome Biology 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.